A SPECTRAL METHOD FOR INTEGRAL FORMULATIONS OF 
MEDIUM-FREQUENCY SCATTERING PROBLEMS 

JOHANNES TAUSCHt 

Abstract. A fast method for the computation of layer potentials that arise in acoustic scattering 
is introduced. The principal idea is to split the singular kernel into a smooth and a local part. The 
potential due to the smooth part is computed efficiently using non-equispaced FFTs, the potential 
due to the local part is expanded as a series in the mollification parameter. The complexity of 
the approach is shown to be 0(n + re 3 log re), where n is the number of degrees of freedom in the 
discretization and re is the wave number. The constant factor in this asymptotic estimate is small 
since no singular surface integrals must be computed. Therefore the method is particularly efficient 
for medium-sized scatterers (50-100 wavelengths) that may have complicated geometry. 
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1. Introduction. It is commonly accepted that the boundary element method 
is an effective approach to solve the Helmholtz equation in the exterior of a scatterer. 
It has the advantage that only the finite boundary surface has to be discretized and 
that the radiation condition is automatically satisfied. In the recent past a variety of 
methods have been developed to handle the dense matrices associated with discretized 
layer potentials. These can be roughly classified into two groups, namely hierarchical 
and grid-based methods. 

Examples of hierarchical methods are the Fast Multipolc Mcthodjn], hierarchi- 
cal matrices and wavelets These methods are based on clustering interactions 
between panels in a hierarchical manner; the larger the separation, the larger the 
clusters. The efficiency and accuracy of these methods depends critically on how 
the cluster interactions can be approximated by low-rank matrices. In the case of 
boundary integral operators associated with the Laplace, Stokes or Lame equations 
asymptotically optimal schemes have been developed. That is, the complexity of a 
matrix- vector multiplication is order n, or order nlog p n, while the convergence rate 
of the discretization scheme is preserved, see, e.g. [T7I[T^1[TU] . 

In the case of the Helmholtz equation, the size of the scatterer, measured in 
wavelengths, is the dominant factor that influences computational cost and accuracy. 
It is well known that in the high-frequency regime large clusters are no longer ap- 
proximated by low-rank matrices, and therefore the efficiency of the aforementioned 
methods breaks down. To overcome this problem it has been proposed to use the Fast 
Multipole Method with diagonal translation operators, \R^. This technique has been 
extended in PITS], 

A different approach that avoids large clusters is to replace the surface distri- 
bution by equivalent sources on a uniform grid. The fast Fourier transform can be 
employed to compute grid potentials efficiently. Since the grid is only accurate when 
the source and the evaluation point are well separated, the nearby interactions must 
be computed directly, by adding up contributions of individual sources. Grid-based 
methods are quite popular even though it appears that the asymptotic complexity 
is generally higher than what can be achieved with hierarchical methods. However, 
in many engineering applications the geometry is complicated and the mesh is rela- 
tively coarse, therefore constant factors often play an important role. Applications of 
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grid-based methods for the Laplace equation can be found in ^] ^] for elasticity in 
[T3*| . For high-frequency scattering this methodology, combined with local high-order 
discretizations have been described in 0|. 

If the scatterer is smooth and isomorphic to a sphere, spectral element methods, 
based on expansion of the solution into spherical harmonics have been shown to be 
successful |E]- 

The approach described in this paper is closer to grid-based methods in that 
FFTs are used to accelerate the matrix-vector product. However, there is no uniform 
grid with equivalent charges. The idea here is to split the Green's function into a 
smooth approximation and a singular, essentially local part. The smooth part of the 
Green's function is replaced by a rapidly converging Fourier series. We will show 
how non-equispaced FFTs can be used to compute layer potentials with such a kernel 
effectively. A similar idea has been applied earlier to the heat equation |10| . 

The local part can be evaluated using expansions with respect to the mollification 
parameter. Thus the computation of the local part amounts to multiplying with a 
diagonal matrix. Since there is no need to compute the nearfield directly, we believe 
that the discussed approach is competitive with hierarchical and grid-based methods. 

We will discuss how the mollification parameter and the number of Fourier modes 
have to be selected as a function of the meshwith and the wave number to obtain 
efficient and accurate schemes. Our analysis is based on bounding the error of the 
bilinear form when the wave number is increased. It should be noted that this does not 
give estimates of the error of the solution. For that, realistic estimates of constants in 
the inf-sup condition are necessary, which are not available. However, we will present 
numerical examples that suggest that the selection strategy of the parameters indeed 
control the error when increasing the wavenumber. 

2. Problem Formulation. For simplicity of exposition, the focus of this paper 
will be on the sound-soft acoustic scattering of an incoming field it lnc of a smooth 
obstacle D C R 3 . The reflected field u is described by the Dirichlet problem to the 
Helmholtz equation with the Sommerfeld radiation condition 



Here, k is the wave number. We assume that the problem is scaled such that the 
scatterer is located inside a cube of side length 1 — d, that is 



where < d <C 1 is a constant. 

A classical approach to treat the Helmholtz problem l|2.1|) is the combined layer 
ansatz of Brakhage and Werner where the scattered field is represented by a 
combination of a single- and double layer potential 



(2.1) 




0, l£l 3 \fl 

xeS:=dD 



(2.2) 



S C [0, 1 - d] 3 



(2.3) 



U 



(x) = {JC - ir]V)cr(x) 7 x£R 3 \D, 



where r\ > is the coupling parameter, a an unknown surface density and 



(2.5) 



(2.4) 
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are the single and double layer operator, respectively. By letting x — > S from the 
exterior of the scatterer, and taking the jump relations of layer potentials under 
consideration, the following boundary integral equation for a can be derived 

(2.6) -a(x) + (K - irjV)a(x) = -u mc (a;), x e S. 

It is well known that l|2.6|) is a well posed problem when rj > 0, see, e.g., [5]. 

To obtain a discretization of l|2.6|) we introduce the space Xh of piecewise poly- 
nomial functions on a triangulation of S. This triangulation is assumed to be quasi- 
uniform and the maximal diameter of a triangle is denoted by h. The nodal basis {<fi}, 
i = 1, ... ,n of Xh consists of functions with local support. Note that n = 0(h~ 2 ). 
The Galerkin approximation is defined to be the function a% S Xh whose residual is 
L 2 (S')-orthogonal to Xh- This leads to the linear system (M + A)x = b, where the 
coefficients of the system matrix and the right hand side are given, respectively, by 

Mij = \ J (Pi(x)<fij(x)dS x 

1 JJ S W y - l V 4n\\x-y\\ nV)<PMdS v dS t , 
ipi{x)u iI,c (x)dS x . 

Since this system is large, iterative methods for its solution must be employed. The 
dominant cost in such a scheme is the multiplication of a vector with the dense 
matrix A, which, if it is done in the obvious way, has 0(n 2 ) complexity. The article 
will discuss a scheme to compute the product approximately with a highly reduced 
flop count. 

3. Splitting of the Helmholtz Kernel. The heart of the method is the split- 
ting of the Helmholtz kernel 

(3.1) G(r) = G s (r) + E s (r) 

into a smooth part Gs(r) and a singular, local part E$(r). Here, <5 is the mollification 
parameter that controls the smoothness of Gg- This splitting results in a splitting of 
the single layer potential 

Vg(x) = $ s (x) +<S> L (x) 

where 



(3-2) &>(x)= / G 5 (x-y)g(y)dS y 

Js 

(3.3) $ L (x) = [ E s (x-y)g(y)dS y 

Js 

The splitting of the double layer potential JCg{x) = ^ s (x) + ^ L (x) is defined analo- 
gously. 

3.1. Smooth Part. The Green's function can be expressed in Fourier space, 

Cq a\ m \ exp(zK||r||) 1 / I ,. . 3 

(3.4) Gir) — — — — I ^ cxpur • u) d w. 

4^||r|| (2tt)3 J w , \\uj\\ 2 -k? 
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Exploiting the spherical symmetry of the Fourier transform leads to 

(3.5) G{r) = l_J^_^_ Mn)dp 

where jo(z) — sin(z)/z is the spherical Bessel function of order zero. The integral in 
(|3.5[) is understood in the sense that the singularity at p = k is circumvented in lower 
complex half-plane, thereby enforcing the Sommerfeld radiation condition, see, e.g., 



The decay rate of the transform at infinity determines the regularity of the kernel. 
The integrand in (|3.5|) is only 0(p~ 1 ) as p — > oo which explains the singularity of the 
Green's function in the origin. A smooth approximation of the kernel can be obtained 
by multiplying the transform with a filter to increase the decay rate at infinity 

(3.6) Gtf ( P ) = _/ H{d(p 2 - K 2 ) )_f!_j (p||r||)dp 

Z7T Jo V / p — K 

where H is the filter. Because of the singularity of the integrand at p = k it is more 
convenient to write the filter in the form as it appears in 1)3.6(1 and not as H(Sp 2 ). 

There are several possible choices for H . If the filter is a rational function, then 
G$ can be expressed in closed form. To that end, write the filter in partial fraction 
decomposition 

H(z) ^ " 



i, z + w k 

k=l K 

where the coefficients Wk and c k are at our disposition. It will become clear later that 
because of the singularity of the integrand it is necessary that 

(3.7) H(0) = 1. 

Basic complex variable arguments show that 

H(z) _ ,X d k 



(3.8) 



where wo = and because of condition (|3.7|l 

(3.9) do = 1 and = — % . 

Substitution of (|3.8|) into l|3.6|l and the change variables p — * ^/~5p leads to 

Gs ^ = ^bn E dk ! ^rhv jo fpfl) dp, 



2^5 to Jo (?-& k JU VV6 

where 

(3.10) k = VSk and Wk 

The integrals in the last expression are of the same form as l|3.5|l . Therefore the 
mollified Green's function has the closed form 

</ 



(3.11) G s (r)= d k 



exp(zK||r||) ^ exp(zw fc ||r||/V^) 



47r||r|| ^ 4tt 1 1 r- 1 
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In the discussion below, it will be convenient to write the decomposition in (|3.11|) in 
the form 

(3.12) G(r)=G((r) + _L E (M 

where E is singular at z = given by 

exp(iwkz) 



(3.13) E(z) = -J2d 



4ttz 

k=l 

A good filter must satisfy two properties. First, it must decay rapidly to ensure 
smoothness, that is, there must be a constant c such that 

(3.14) \H(z)\ < cmin(l,z-«) . 

Second, the kernel E must decay exponentially away from the origin, this is why it 
will be referred to as the local part. To ensure the latter condition, it is necessary 
that Im(wk) is bounded away from zero as k — > oo. If Wk is real, then l|3.10|l implies 
that 



(3.15) y/5 < - 



mm \Wk\- 

K k>0 

The latter condition implies that k is bounded as K — > oo. 
An example of a filter that satisfies (|3.14|) is given by 

(3.16) ^)=n^b 

This filter has poles Wk = iV~k, k = l,...,q. Condition l|3.15|l is equivalent to 
y/S<l/K. The filter 

(3.17) H{z ' 



(! + «)« 



also satisfies the decay property. Because of the repeated poles, the smooth part 
corresponding to this filter is slightly different from l|3.11|l : 

, , exp(i«||r||) exp(iw||r||/v / ^) 
Gs{r) — 



**\\r\\ \VS 

Here p(-) is a polynomial of degree q — 1 and w — yl — Sn 2 . As with the previous 
filter, condition 1|3.15J) is VS < 1/k. 

3.2. Local Part. The smooth part is a good approximation of the actual Green's 
function if S is small and r is large. In the neighborhood of the origin the two functions 
are very different and therefore the contribution of the local part must be accounted 
for. In this section we show that the local part has an expansion with respect to the 
mollification parameter VS and show how to compute the expansion coefficients. 

Since condition l|3.15(l implies that the coefficients Wk in Ij3.13|l have a positive 
imaginary part, the function Eg decays exponentially away from the origin. We 
introduce the smooth cut-off function \ for some < v < 1 which is small enough 
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such that the surface has a parameterization of the form y(t) = x + At + nh(t) in the 
^-neighborhood of x. Here n is the normal of the surface at the point x, A G E 3x2 
has two orthogonal columns that span the tangent plane at x and h(t) = (3(|i| 2 ) is 
some scalar function in t £ M 2 . The local single-layer potential $s(x) hi l|3.3[) can be 
written in the form 



$ S (x) = I E s {x - y)g{y) dS y 

Es(x - y)Xu(x - y)g{y) dS y + O (exp (-^)) 

(3.18) = J ^ E s (t)g(t)d 2 t + O (exp . 

Here, E${t) = E$(x — y(t)), g(t) — Xu(x — y(t))g(t)J(t) and J(t) is the Jacobian of the 
parameterization. For simplicity of the argument we assume that the function h(t) in 
the parameterization of the surface is analytic, that is, 

(3.19) h(t) = M"- 

\a\>2 

Thus there are there are C°°-functions H n such that 

oc 

(3.20) r(t):=||x-y(t)|| = ||i||X;ilC-ffn(t) 

n=0 

where t := t/\t\ and H (t) = 1 and Hi(i) = 0. Substituting IpHH)) into (|3~T%|) leads to 

(3.21) ^) = ^/ R ^(^)^)^ + 0(exp(-9) 

(3.22) =VsJ E ^t\^{VS\t\) n H n (^j g(V5t) d 2 t + O (cx V 

where the second integral is the result of the change of variables t i— ► The 
integral as a function of VS is C°°, and can therefore be expanded in a Taylor series. 
The expansion coefficients are derivatives of the integral with respect to Vo. We see 
that 

(3.23) $ s (x) = Si$ g(x)+O(5i) 



where 



*»=/ E(\\t\\)dh= ^±^. 



k=l 

The double layer potential is given by 

The second factor of the kernel can be expanded in a similar manner as l|3.2U|l . we 
find that 

^ ', % = -lltll (/i 02 cos 2 6> + /iiicos(9sm(9 + /i 20 sm 2 0) + 0(\\t\\ 2 ) 
\\x-y\\ 
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where 9 is the angular coordinate of t and the h's are from expansion l|3.19fl . Pro- 
ceeding in a similar manner as for the single layer operator the following expansion 
follows 

(3.24) 9s{x)=8*y g(x)+O(8i) 
where 



fc=i 



#o = --(ho2 + h 20 V — 



4. A fast algorithm for smooth, periodic kernels. We describe a fast al- 
gorithm for the smooth part of the single layer in (|3.2|) . which is based on Fourier 
analysis. Modifications for the double layer are minimal and mentioned at the end 
of the section. For rapid convergence of the Fourier series it is necessary to multiply 
the smooth kernel with a sufficiently smooth cut-off function G s := \Gs that is unity 
inside the cube [— 1 + d, 1 — d] 3, and vanishes outside [— 1, l] 3 . Recall that we assumed 
in 12.2f> that the surface is contained in [0, 1 — d] 3 , thus the cut-off function has no 
effect in the integral, and the smooth part is given by 

(4.1) <P s (x) = [ G s (x-y)g(y)dS y: x e S. 



J.S 

The kernel G s can be approximated by the truncated Fourier series Gn 

(4.2) G N (r):= ]T G k exp(7rz k T r) , r E [-1, l] 3 

||fc|| <JV 

II II oo — 

where the summation index A; is in Z 3 . The resulting approximate potential is given 

by 

(4.3) $n(x) = / G N (x - y)g(y)dS y = V] exp(wik ■ x)d k 
where d k = Gkgk and 



11*11 ao^ 



(4.4) G k = l [ exp{-nik-r)G s {r)d 3 r, 

8 -A-M] 3 

(4.5) dk= I exp(-7ri/c • y)g{y)dS y . 

Js 

In case the smooth part of the double layer is to be calculated, the coefficients 
must be replaced by 



(4-6) % 



f d 

J s dn~ exp ( _7rifc ■ y)9(y)ds y . 



To simplify the discussion, our notations will not distinguish between the coefficients 
in l|4.5|l and (|4.6() . In summary, the potential computation of the potential due to the 
smooth parts consists of three stages, in (14.311 . 

1. Compute the Fourier coefficients 'gk in l|4.5jl . 

2. Multiply d k := G k g k for Hfc^ < N. 
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3. Evaluate the Fourier series (|4.3|) for x G S. 
The choice of the truncation parameter N depends on the wave number and the 
mollification parameter and can be much smaller than the linear system size n. In 
Section|Slthe exact dependence will be investigated. Stage 2 obviously involves 0(N 3 ) 
operations, the other two stages can be executed efficiently using non equispaced Fast 
Fourier Transforms. This will be discussed next. 

4.1. Computation of the t/fc's. In this section we describe how FFTs can be 
used to efficiently compute the Fourier coefficients of the function g. To that end, the 
three-space is divided into small cubes C/, I — (h,h,h) & Z 3 ,0 < lj < N. These 
cubes have centers xi = l/N and side length 1/N. Note that ./V is the same as in ()4.2[1 
and therefore the cubes get smaller if more terms in the Fourier series expansion of 
the Green's function are retained. Because of assumption l|2.2|> S is contained in the 
union of all cubes and set Si = Ci n S to denote the piece of the surface that intersects 
with the Ith cube, c.f. Figure |4~T1 




(4.7) 



Fig. 4.1. Two dimensional illustration of the geometry. 

From (|4.5|l it follows that the Fourier coefficients of g can be written as 



9k = ^2 exp 

IUII_<JV 



N 



exp(-7rifc • (y - xi))g{y) dS y 



s, 



The frequency and the spatial variable in the integral can be separated using the 
Jacobi- Anger expansion 



exp(-i£t) = J2(-i) v (2p + l)j„(£)P v (t), -1 < t < 1, 



!/=0 



see, e.g., ^21- Here, j v (-) is the spherical Bessel function of order v and P v (-) is the 
Legendre polynomial of degree v. This formula generalizes to the three-variate case 
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and can be applied to the integrand in 14.7)1 

(4.8) exp(-7rifc • (y - a;,)) * E H) H (2a + l)j a (nkH)P a (^j^) 

\a\<p 

where p is the expansion order, H = 1/(2N), a — [a±, 012, 0.3) is a multi-index, 
\a\ — ai + a.2 + «3, j a {x) — j ai {xi)ja 2 ( x 2)ja 3 {x3) and P a {x) is defined similarly. 
Substitution of 1)4. 8fl into (|4.7I) leads to the approximation 



where 

(4.9) mf(g) = J ^ P a { v -^L) g ( y ) dS y 

is a moment for which exact formulas can be derived if the function and the surface 
are discretized. In particular, if g is a piecewise polynomial, then the moments are 
linear transformations of the coefficients of g corresponding to the nodal basis. The 
matrix that maps the coefficients to the a-th moments is denoted by M a . The number 
of nonzero entries in M a is n. 

In matrix form, the (approximate) coefficient vector g is given by 

(4.10) g= K a FM a g, 

\a\<p 

where F is the 27V-long three-dimensional discrete Fourier transform, g the vector of 
coefficients of g and K a is a diagonal matrix with the factors (— i)' Q ' (2a + l)j a {-nkH). 
The computation of g involves (p + l)(p + 2)(p + 3) /6 FFTs. In Section[S]we will show 
that it suffices to use a small value of p. 

4.2. Evaluation the Fourier series. In the Galerkin discretization, the i-th 
component of the matrix-vector product <I>i is the inner product of the potential <& 
in 1)4. ljl with the i-th nodal basis function ipi. For the fast method, the potential 
is replaced with the approximated potential $>n in 1)4.3)1 . In order to evaluate the 
potential efficiently, the Jacobi-Anger approximation ()4.8)) is used again, in a very 
similar manner as in the previous section. This is shown in the following computation 



$ 4 = / ipi(x)<$> N (x)dS x 
s 

E exp ( *i 1 j f exp(irik ■ (y - xi))<pi(x)dS x dk 



~E E exp(^^ji^(2a+l)j a (TTkH) m f(cpi)d k . 
In matrix notation, the above can be written as 

(4.11) $ = E mIf*kJ. 

\a\<p 

Hence (p+l)(p+2)(p+3)/6 FFTs are necessary to compute the vector Furthermore, 
it is evident that the operation ()4.11)) is the adjoint of operation ()4.1U)) . 
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5. Error Analysis. In this section we derive estimates for the error of the bi- 
linear form introduced when the mollified kernel corresponding to the single layer 
operator is replaced by the truncated Fourier series expansion. To keep the technical 
level of the discussion at a minimum we omit the discussion of the double layer oper- 
ator, because it is completely analogous to the single layer. The main concern is the 
situation where k — > oo and our goal is to determine N and p as a function of the 
wavenumber such that the resulting error remains bounded. 

It is straightforward to see that 

(/, (A - A N )g) = dk ^9k 

l|fe|lco>W 

where A is the surface integral operator with kernel G s , An its Fourier series approx- 
imation, Gk are the Fourier coefficients of G and gk are the Fourier coefficients of 
surface distributions as defined in 1)4. 5[l . The obvious way to estimate the error is 

\(f, (A-A N )g)\ < Yl \ d * 

llfcll >JV 

(5-1) < e \ d * 

It is not possible to work with ^-estimates of because / can be regarded as 
surface-delta function in R 3 which cannot be bounded by the L 2 -norm. Thus the 
Fourier coefficients of the kernel must be estimated in the /i-norm which amounts to 
an estimate in the L^j -^-norm. 

The Fourier coefficients can be related to the derivatives of the function with the 
standard integration-by-parts argument. Since the kernel is a spherically symmetric, 
three-variate function, it is convenient to work with the Laplacian and the Green's 
formula. Because of Aexp(i/c • r) = — \\k\\ 2 exp(ifc • r) it follows that 

Gk = \l G s (r) X (r)d 3 r 
° J[-i,i] 3 

= 1^1 A™(G,(r) X (r))d 3 r = -^ r [A^ X )] fe , 
8||fc|| J[-i,i]i ||fc|| 

for any integer m for which the right-hand side is defined. With this estimate at hand, 
one obtains 

E N<(e|[^x)l| 2 V ( E -p) 

||fe|| oc >A r \feez 3 / \||fe|| OQ >A r H^ll / 

(5.2) <cNl~ 2m \\A m (Gsx)\\ L z 

where the first step follows from the Cauchy-Schwarz inequality and the second step 
follows from Parseval's equation. Since the mollified Green's function gets more 
peaked in the origin as S — > the norm of Gs cannot be treated as a constant. 
Therefore N must be linked to 6. Unfortunately, the right-hand side in (|5.2J) involves 
the product of the Green's function with the cut-off function. The product rule leads 
to estimates that involve factors which depend on m and which are difficult to control. 



sup 

u 



SUP \dk\ 

I. 



\\f\\L2(S)\\9\\L 2 {S)- 
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Furthermore, the argument in (|5.2|l assumes that the cut-off function has the same 
regularity as Gg, which, as numerical experiments suggest, is not necessary. 

The following discussion presents a refined analysis intended to obtain more re- 
alistic error estimates. 

5.1. Estimates of the Derivatives. Derivatives of the mollified Green's func- 
tion can be obtained either from the Fourier transform (|3.6|1 or the closed form Q3.11|l . 
The Fourier integral leads to estimates that display the dependence on the order of 
differentiation more clearly. For the subsequent error analysis it suffices to work with 
powers of the Laplacian. Since the Bessel function satisfies the Helmholtz equation it 
follows from (|3.6(l that 



for any integer m < q. Because of the assumption (|3.14|) , the filter decays only if the 
argument is larger than unity. This condition leads to 



(5.3) 




(5.4) 




The definition of po immediately implies that 



(5.5) 



Pa = —r= V 7 ! + K 2 , 



i.e., po < en. 

Lemma 5.1. The following estimates hold 



(5.7) 



(5.6) 




Proof. Break integral (|5.3() into three parts 



A m G s (r) 



(-1) 



m 



(Ji + h + h) 



where 




and 



f{p) := p 2m H 



12 



J. TAUSCH 



Since Sp 2 , < c the derivative of f(p) is bounded by \f'(p)\ < cmp^" 1 1 for < p < p a . 
Thus the first integral can be estimated as follows 



\h\< max |/'(p)| r -B^\j {p\\r\\)\dp 

0<P<Pa J p + K 



p 2 ™- 1 [»» p 



< cm— / dp 

r Jo P + K 

2m 

(5.8) < cm^. 

\\r\\ 

The integral in 1% can be computed in closed form. The resulting expression involves 
the integral sine and cosine functions and are easily shown to be uniformly bounded, 
thus 



(5.9) \h\<c- 



r 



The integration in I3 is over the interval where the filter is decreasing. Therefore 



\h\ < TT-TTT^ / 7^E .9 ~ .9 d P 



r 



1 

< 11 11 ; max 
~ \\r\\5i po<p 





1 


(p 2 — K 2 )? p 2 — K 


p2m+l 


poo 


(p 2 - K 2 )1 


J pa 



Po 

1 



p 2 — K 2 



dp. 



The function to be maximized is monotonically decreasing; the integral can be com- 
puted in closed form and estimated by c/k. Thus ^3 can be estimated by 

(5.10) |/ 3 | < -^Pl m+1 . 

n\\r\\ 

Estimate (|5.6|) is immediate from (|5.5|) . I|5.8|l . (|5.9I) and (|5.1U|) . Estimate l|5.7fl follows 
from 1|5.6[) because the 1/r singularity cancels upon integration. □ 

Using very similar arguments as in the previous proof, the first and second deriva- 
tives of powers of the Laplacian can be estimated. We only state the result. 

Lemma 5.2. For \a\ < 2 we have 



(5.11) \\d a A™G s (r)\\ L2 <c—- T (l + k 



, 171 ft , 

L f-l,l]3 - ^S m+ ^ 



5.2. Approximation Analysis of the Fourier Series. Our goal is an estimate 
the Fourier truncation error in the spirit of i|5.2f> . that does not involve high-order 
derivatives of the cut-off function. 

Lemma 5.3. If x € C 3 then the Fourier coefficients of the kernel are given by 

Gk = ([^Gs] k + 2m[WA"^G 5 ■ V X ] fc 

rn — 1 1 _ 

1 = \\ k \\ l = |e{1.2} 

|/3|=4-|a| 

for any integer m for which the right-hand side is defined. The coefficients satisfy 
<4 )/3 G{0,l,2Z,4Z}. 
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Note that the above expression only contains derivatives of the cut-off function up 
to order three. Furthermore, the remainder (i.e., the sum over I), contains derivatives 
of Gs which are at least two orders lower than the power of l/||fc||. This will be 
essential for the subsequent error analysis. 

Proof. A simple application of the product rule shows that 

A(G 5X ) = AG SX + 2VG 5 • V x + G S A X , 
A(VG 5 ■ Vx) = VAG 5 ■ A X + 2tr(G' s 'x") + V X ■ VAGa. 

From integration by parts it follows that 

Q k = -^[MGT X )] k = -jp ([AG^] fe + 2[VGTvx] fc + [G^lJ ■ 
||fc|| \\k\\ V / 

Repeating this argument for the first two terms and leaving the third term unchanged 
leads to 

G k = -lj ([A^G7x] fe + 4[VAgTvx], 
ll fc ll v 

+ T^P (t^A^fe + + 2[VG7^VA x ] fc ) + -^[OAx} k . 

\\k\\ v / \\k\\ 

By induction one finds that 

& = ^pr ([A^G7x] fc + 2m[VA"^ • V X ] fe ) 

m — 1 1 _ 

+ E ^1+2 {i^K x ] k + ii[ti-A^G" X '% + 2/[va'-^gT: va x ; 

1=0 11*11 

which is the assertion. □ 

Combining this result with the estimates of the derivatives of Gs leads to the next 
theorem 

Theorem 5.4. If \ € G 3 , /, g e L 2 {S) and N is chosen such that 
(5.12) A:=(i±|!) 4 <l, 

then for any integer < m < q the approximation error of the truncated Fourier 
series is bounded by 

|(/, (A A N )g)\ < c (m 2 Ni\ m + AT*) \\f\\ LHS) \\g\\ LHsy 



This theorem suggests how N must be selected to control the error as k — * oo. 
Recall that condition (|3.15(l implies that 6 ~ 1/n 2 to ensure that k is bounded. 
Because of Q5.121 one has to select N such that N ~ 1/ V~5- Since q is free, the product 
can always be controlled by letting q increase as N increases. The influence of 
q on the computational cost is negligible, therefore the error can be controlled with 
complexity 0(n + k 3 log k). 



14 



J. TAUSCH 



Proof. Using the previous lemma it follows that 

1 



E l^k < E _>,„ 

||fe|| DO >iV l|fc|loo> Ar " " 

m— 1 



[A m G 5 ] k + 2m[VA m ^Gi • V x ] k 



E E E 

/=0 l«|e{i,2} llfcll >at 

|/3|=4-| a | "°° 



1 



2i+2 l"a,/3 



[d"A. l -iGsdPx\k 



Using Cauchy-Schwarz and Parseval in a similar manner that lead to estimate (|5.2|l 

1 



E & 

l|fc|loo> JV 



< — — ^ (||A m G A - X || + IIVA™- 1 ^ • Vx||) 

rn—1 

+ E^itt E Kpiiifl-A'-^^xii. 

i=0 " | Q |e{i,2} 

|/3|=4-|a 

The derivatives in the above expression can be estimated using the inequalities (|5.7J) 
and H5.11f> derived in the previous section. Recalling the definition of A in (|5.12|) , this 
leads to 

rn 1 

\Gk\< cNim 2 \ m + cN-iJ2 Pxl - 

l|fc|| 00 >A r '=0 

Since A < 1 the sum is bounded independently of m. Combining the last inequality 
with 1)5. 1[) completes the proof. □ 

6. Error Analysis of the non-equispaced FFT algorithm. The error anal- 
ysis of the previous section is not complete since in the non-equispaced FFT algorithm 
the complex exponential function is approximated by the truncated Jacobi- Anger ex- 
pansion. It is therefore important to know how the highest retained order p must be 
selected as a function of the wave number. This will be determined in this section. 

6.1. Error of the multivariate Jacobi- Anger approximation. The Jacobi- 
Anger expansion is an expansion in Legendre polynomials. Therefore the error of the 
multi-variable truncated expansion in (|4.8(l is given by 



(6.1) 



e£(i):=e fc (t)-eS(t)= ]T £ e a L a (t) 

n=p+l |a|=n 



where e k (t) = exp(iH k-t), e k denotes the truncated expansion and the coefficient e c 
has the form 



-i,ip 



e k (t)L a (t)d 3 t. 



(-1)" 
a\2 a 



[-1,1] 



{l-t 2 ) a d a e k {t) d A t. 



The latter form follows from the Rodrigues formula and integration by parts. It is 
useful for estimating the magnitude of the coefficient 



|e a | < 



1 



a\2 a 



[-1,1]' 



(l-t i ) a d*t\\d a e k \\ 00 < 



1 



a!2 Q 



(vriJfc)". 
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Since |i n (*)| < 1 f° r |*| < 1; the truncation error of l|6.1|) can be bounded as follows 



%®\ = E E w» 

n=p+l 



* E 

n— p+1 



n— p+1 



E 

|a| — n 

1 1*11? 



The last step is an application of the multivariate binomial formula. Since the last 
equation is the remainder of the Taylor expansion of the exponential function, we 
have the bound 



(6.2) 



! / 7T 



< 



If in (|6.1[) the Taylor series instead of the Jacobi-Anger is employed, then a very 
similar analysis shows that 



P+1 



Thus the Jacobi-Anger expansion is significantly more accurate for large values of /c. 

6.2. Error of the non-equispaced FFT. In the non-equispaced FFT algo- 
rithm, the kernel Gat in 14.2(1 is replaced by the kernel 



(6.3) G p N (x, x') = Yl ^ ex P 



k-(l- V) 
N 



ek{ti)e k {t' v ), x e Ci,x' e C v 



HfclU<JV 

where ti — (x — xi)/H, t\, = (x! — xi>)/H. Thus the error is given by 
G N {x - x) - G p N (x, x') = 



E Gk exp 



fc • (Z - V) 

' N 



(ek(ti)ek(t'i,) + ek(U)ek(t' v ) + ek(U)e k {t' l ,) 



From l|6.2|l and ||fc||-, < v3||fc|| the estimate 



(6.4) G N (x-x')-G p N (x,x') 
follows. We use the identity 



< 



(p + 1)1 1 ~W 



p+i 



E | e *n*ir +1 



G k 



[A 2 (Gx)] A 
llfcll 4 



and recall that from Section [S] it follows that 

l|A 2 (G«)ll L2 < c <r 2 . 



16 



J. TAUSCH 



Continuing with estimate 1|6.4H gives 
G N (x-x')-G p N (x,x')\ 

v 7 V / \ll*IL<w / \llfclUSJv / 

We have almost completed the proof of the following theorem 

Theorem 6.1. Let A P N the integral operator that has kernel G P N then 

\(f,(A N -A» N )g)\ < (^j N-i \\f\\ LHS) \\g\\LHsy 

The most important conclusion from this result is that the order p in the Jacobi- 
Anger approximation does not have to be increased as N — > oo. 

Proof. Elementary integral calculus and estimate (|6.5|l together with the fact that 
S ~ iV -2 imply that 

|(/, (An -A p N )g)\ < max \G N (x - x') - G p N (x,x')\ ||/|| l i (S )II.9|Ili(s) 

□ 

7. Numerical Examples. We have implemented the method to verify the the- 
oretical estimates. In this implementation the Fourier coefficients are computed nu- 
merically. The method used for this task is completely analogous to the computation 
of the Fourier coefficients of a surface density described in Section 14.11 The only 
difference is that the moments in (|4.9|) are replaced by the moments of the Green's 
function 

mt(G) = jf P a {U^l)G s {y) X {y)d?y 

where C; = xi + H[— 1,1] . These moments are computed using Gauss quadrature. 
The analysis of the error introduced by computing Gk numerically parallels the dis- 
cussion of Section l6~2l and is therefore omitted. 

In the first example we compute the farfield pattern when the unit sphere is hit 
with a plane wave. This is done by solving integral equation Q2.6|l with piecewise 
constant elements combined with the spectral method. The farfield is computed from 
the density using the formula 

a(x) = i / exp(—inx • y)(r] + kx ■ n y )a(y) dS y , ieS, 
Js 

see, e.g., ^21- Because of the spherical symmetry the solution a as well as the farfield 
can be expressed in closed form. The coupling parameter in (|2.3|) is rj — k/2 and the 
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linear system is solved with GMRES without any preconditioning. Tables [7TI and I7T1 
display the relative errors of the L 2 - n orm of a(x) when increasing the size of the 
sphere measured in wavelengths. The results show that the error remains bounded 
(actually decreases somewhat) when N ~ n and 5 ~ 1/n 2 and are therefore in good 
agreement with the theoretical estimates. The truncation parameter of the Jacobi- 
Anger expansion in (|6.1[1 is always set to p = 4. In these experiments the meshwidth 
is proportional to the wavelength, which is reflected in the fact that the number of 
panels n is quadrupled in every line. We have implemented both filters (|3.16[1 and 
l|3.17[) and set q — 5. The displayed results are for (|3.17l) . but the results for the other 
filter are only marginally different. The timings displayed are the time per iteration 
and the total time, which also includes the time to compute the Fourier coefficients 
Gk- The cpu is a 3.6 ghz Intel Xeon processor. The time per iteration increase by a 
factor somewhat larger than eight when doubling the wavenumber, which agrees well 
with the k 3 log k complexity estimate. The code stores the Fourier coefficients Gk and 
gk, the moments m a (g), and as well as the orthogonal basis of the Krylov subspace 
generated by GMRES. For the size of problems computed, the basis consumes the 
largest portion of the overall memory usage. Since this part grows roughly like k 2 , 
the growth rate of the overall storage appears slower in the Table than the asymptotic 
k 3 log K estimate. 



n 


N 


size 


8 




its 


mem 


time/itr 


time 


error 






(A) 








(MB) 


(sec) 


(sec) 




5120 


16 


6.25 


1.00x10" 


-4 


11 


8.0 


0.4 


5 


0.076 


20480 


32 


12.5 


2.50x10" 


-5 


12 


32.8 


2.1 


32 


0.044 


81920 


64 


25 


6.25x10" 


-6 


15 


139.8 


18.3 


337 


0.039 


327680 


128 


50 


1.56x10" 


-6 


18 


623.9 


159.2 


3355 


0.034 


1310720 


256 


100 


3.91x10" 


-7 


22 


2981.7 


1386 


34414 


0.031 



Table 7.1 
Results for the sphere. Lower accuracy. 



n 


N 


size 


6 




its 


mem 


time/itr 


time 


error 






(A) 








(MB) 


(sec) 


(sec) 




5120 


16 


3.13 


1.00x10" 


-4 


8 


8.0 


0.4 


5 


0.044 


20480 


32 


6.25 


2.50x10" 


-5 


11 


32.8 


2.2 


32 


0.021 


81920 


64 


12.5 


6.25x10" 


-6 


12 


139.8 


17.9 


276 


0.011 


327680 


128 


25 


1.56x10" 


-6 


15 


623.9 


159.3 


2880 


0.0071 


1310720 


256 


50 


3.91x10" 


-7 


18 


2981.7 


1378 


28714 


0.0051 



Table 7.2 
Results for the sphere. Higher accuracy. 



To illustrate that the technique discussed in this paper can be used for very general 
scatterers we include the Boeing 747 example shown in Figure 17.11 The surface of 
the airplane is assumed to be sound soft. The geometry is given by a list of vertices 
and triangular panels which can be downloaded from the internet. There are 556552 
panels, and further information, such as parameterizations, are known. We ignore the 
fact that there are edges and conical vertices in the geometry and set the curvature 
term in (|3.24() to zero. 

Figure 17.21 compares the density for N = 128 and N — 256 Fourier modes. 
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Fig. 7.1. 3D rendering of the airplane. 

Since it is hard to spot differences in the two solutions, it appears that already the 
smaller value of N will give an acceptable accuracy in many applications. The size of 
the scatterer in this problem is about 45 wavelengths, the memory allocation of the 
smaller problem is 906MB and the cpu time is 4355 seconds. 

Figure 17751 displays the solution for 90 wavelengths and N — 256. The memory 
allocation is 1583 MB and the cpu time is 36361 seconds. 




Fig. 7.2. Comparison of the density (imaginary part) for N = 128 (left) with N = 256 (right); 
45 wavelengths 

8. Conclusions. We have presented a method for the computation of scattered 
fields that has 0(n 3 log k) complexity when the meshwidth is proportional to the wave- 
length. Since n ~ k 2 the asymptotic estimate is not optimal, but because of small 
constants we have been able to solve lOOA-problems in eight to nine hours. Most of 
the cpu time is spent evaluating the sums in (|4.10|) and l|4.11|l . Since this part is em- 
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Fig. 7.3. The density (imaginary part) for 90 wavelengths 



barrassingly parallel one can expect almost optimal speed up on distributed memory 
multiprocessor machines. The approach generalizes to electromagnetic scattering. 

9. Acknowledgement. The author obtained the panel description file of the 
airplane from the website www.3dcafe.com. 
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